Oscillatory Notch pathway activity in a delay model of neuronal differentiation 
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Lateral inhibition resulting from a double-negative feedback loop underlies the assignment of 
different fates to cells in many developmental processes. Previous studies have shown that the 
presence of time delays in models of lateral inhibition can result in significant oscillatory transients 
before patterned steady states are reached. We study the impact of local feedback loops in a model 
of lateral inhibition based on the Notch signalling pathway, elucidating the roles of intracellular and 
intercellular delays in controlling the overall system behaviour. The model exhibits both in-phase 
and out-of-phase oscillatory modes, and oscillation death. Interactions between oscillatory modes 
can generate complex behaviours such as intermittent oscillations. Our results provide a framework 
for exploring the recent observation of transient Notch pathway oscillations during fate assignment 
in vertebrate neurogenesis. 
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I. INTRODUCTION 



As in many electronic circuits, classes of oscillators and 
switches are fundamental elements in many gene regula- 
tory networks [l^, H^. In particular, a double- negative 
feedback loop comprising two mutually repressive com- 
ponents is known to be capable of functioning as a toggle 
switch, allowing a system to adopt one of two possible 
steady states (corresponding to cell fates) [1, [l] . In the 
context of developmental biology, such bistable switch 
networks can operate between cells, and are believed 
to drive cell differentiation in a wide range of contexts. 
However, in naturally evolved (rather than engineered) 
gene regulatory networks, double- negative feedback loops 
rarely exist in a "pure" form, and interactions between 
loop components and other network components often 
result in sets of interconnected feedback loops. Further- 
more, if loop interactions involve the regulation of gene 
expression, then interactions are delayed rather than in- 
stantaneous. The present study investigates the dynamic 
behaviour of a double-negative feedback loop when the 
nodes of the loop are involved in additional feedback 
loops, and when the regulatory steps constituting the 
resulting network involve significant time delays. 

A particularly well documented example of a biological 
double-negative feedback loop is centred on transmem- 
brane receptors of the Notch family. Notch signalling, 
resulting from direct interaction with transmembrane lig- 
ands of the DSL (Delta, Serrate and Lag-2) family on 
neighbouring cells mediates an evolutionarily-conserved 



lateral inhibition mechanism that operates to specify dif- 
ferential cell fates during many developmental processes 

0, [n, El, m. Although gene nomenclature varies be- 
tween different organisms, a core circuitry — the neuro- 
genic network — underlying lateral inhibition can be iden- 
tified, and is illustrated schematically in Fig. [I][l,[l3|. In 
brief, signalling between neighbouring cells is mediated 
by direct (juxtacrine) interactions between Notch recep- 
tors and DSL ligands. A double-negative feedback loop 
is formed by repression of DSL ligand activity by Notch 
signalling in the same cell (cell autonomous repression) — 
Fig. [TJa). Mathematical models of such a spatially- 
distributed double-negative feedback loop are capable of 
generating robust spatial patterns of Notch signalling in 
populations of cells [1| . 
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FIG. 1: The neurogenic network models in a two-cell repre- 
sentation. The letters on the nodes correspond to the key 
classes of protein involved in the network — D: DSL ligand, 
N: Notch receptor, H: Hes/Her proteins. A: proneural pro- 
teins (e.g. Achaete, Scute, Neurogenin). Edges represent two 
types of directional interactions: activation (— >) and repres- 
sion (— o). (a) The pure double-negative feedback loop involv- 
ing only DSL ligand and Notch receptor [2^ . (b) A more 
detailed model that incorporates Hes/Her negative feedback 
and proneural positive feedback [TtI . 
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In many developmental settings, the level of Notch sig- 
nalling regulates the fate adopted by a cell by acting as 
an input to a cell autonomous bistable switch formed by 
one or more proneural genes (such as achaete and scute in 
Drosophila and neurogenin and atonal in vertebrates) 
The basic principle underlying this switch is the ability of 
the protein products of the proneural genes to positively 
regulate transcription of proneural genes, resulting in a 
direct positive feedback loop. In many systems, including 
the developing nervous system. Notch signalling regulates 
the proneural switch via regulation of the expression of 
proteins of the Hes/Her family. These proteins act as 
transcriptional repressors, and can repress their own ex- 
pression and interfere with proneural self-activation 
Furthermore, proneural proteins can also positively reg- 
ulate the expression of DSL proteins, forming a complete 
circuit of interactions as shown in Fig. [Ijb) . Considered 
as an intercellular signalling network — which we shall re- 
fer to as the neurogenic network — this circuit comprises a 
spatially-distributed double-negative feedback loop with 
additional local positive and negative feedback loops. 

A detailed mathematical model of the neurogenic 
network, incorporating Hes/Hcr negative feedback and 
proneural positive feedback, has been studied by Meir 
et al. [itI , who showed by computer simulation that the 
network is capable of generating spatial patterns of Notch 
signalling in populations of cells. The models of Collier 
et al. Q and Meir et al. incorporate the implicit as- 
sumption that all interactions are non-delayed. However, 
in reality the basic production mechanisms that regulate 
gene expression (gene transcription and mRNA transla- 
tion) are associated with time delays [3, [l^ . Incorpo- 
ration of explicit time delays in the pure double-negative 
feedback loop shown in Fig. [Ija) results in competi- 
tion between dynamic modes, with stable spatial pattern- 
ing typically preceded by significant oscillatory transients 
[2a |. In a biological context, such transients would play 
an important role in delaying the onset of cell differen- 
tiation in a developing tissue. Delays can also generate 
oscillatory dynamics in models of Hes/Her negative feed- 
back loops d, [ll|, [H, [l^ , and such oscillations have been 
observed experimentally 0, H, [I!] . 

As predicted by mathematical models [ll|, H^, oscil- 
latory expression of DSL ligands, Hes/Her proteins and 
proneural proteins has recently been observed in neural 
precursors in the developing mouse brain (23 |. Further- 
more, these oscillations have been predicted to play a 
central biological role in delaying the onset of neural dif- 
ferentiation [2J]. In principle, these oscillations could be 
driven by the cell autonomous Hes/Hcr negative feedback 
loop, with Notch signalling providing coupling between 
cells jnl, or by the double- negative feedback loop cen- 
tred on the DSL-Notch interaction [2g]. In the follow- 
ing, we investigate the interplay of local and intercellu- 
lar feedback loops in models of the neurogenic network, 
using a combination of linear stability analysis and nu- 
merical simulations, emphasising the dynamical effects of 
the multiple time delays in the network. We study prin- 



cipally the case of two coupled cells, since this captures 
the essential features of oscillator synchronisation and 
cell state differentiation. We also show how our results 
extend to larger populations of coupled cells. 

II. THE FULL HES/HER-PRONEURAL 
MODEL AND ITS DISSECTION 

In Fig. [Ijb), positive regulation of Hes/Her (H) by 
proneural protein (A) in the adjacent cell, mediated by 
DSL-Notch signaling, can be considered simply 
cade of three low-pass filters. This simplification allows 
the model in Fig. [Ijb) to be reduced to the model in 
Fig. [l{a), where r denotes time delay, and / and g rep- 
resent generic increasing and decreasing functions, re- 
spectively. In this model, referred to hereafter as the 
full Hes /Her-Proneural model, Hes /Her proteins repress 
proneural proteins in the same cell (1 or 2), while proneu- 
ral proteins activate Hes/Her proteins in the adjacent 
cell. This main loop is supplemented by the two lo- 
cal loops: Hes/Hcr auto- repression and proneural auto- 
activation. Each interaction is not instantaneous but in- 
volves a time delay, typically of the order of minutes to 
tens of minutes [1, [ill, [Hj Ell • The delays in the cell- 
autonomous regulatory steps {T2-Ti) originate predomi- 
nantly from processes associated with gene transcription, 
whereby the DNA sequence of a gene is transcribed into 
a corresponding mRNA molecule, while the delay in the 
non-cell-autonomous interaction (ri) represents in addi- 
tion processes involved in DSL-Notch signalling, such as 
protein processing [2l| . 



(a) (b) g. 




(d) g.ts (e) 
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FIG. 2: The neurogenic network models examined in the 
present study, r: time delays, /: a general increasing func- 
tion; g, h: general decreasing functions. With an appropriate 
choice of time delay (ti), the model in (a) provides a simpli- 
fied representation of the full model shown in Fig. [TJb) . The 
models in (a)-(d) have the potential to exhibit both oscilla- 
tions and diff'erentiation, while the model in (e) exhibits only 
differentiation. 

The models in Fig.[2jb)-(e) are obtained by sequential 
reduction of the full neurogenic network. When Hes/Her 
feedback is negligible, the full Hes/Her-Proneural model 
in (a) becomes the model in (b); while when proneu- 
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ral feedback is negligible, the model in (a) becomes the 
model in (c). When both local loops are negligible, mod- 
els (b) and (c) can be reduced to the model in (d). Fi- 
nally, because two sequential repressions function as a 
net activation, the model in (d) can be further reduced 
to the model in (e). Such simple network elements ap- 
pear repeatedly in gene regulatory networks (and possi- 
bly in other biological and non-biological networks), and 
are examples of what are often called "motifs" The 
models in Fig. [D^b)~(e) are called hereafter: (b) the auto- 
activation two-Proneural model; (c) the auto-repression 
two-Hes/Her model; (d) the non- autonomous Proneural 
model; (e) the auto-activation single Proneural model. 



III. MATHEMATICAL REPRESENTATION OF 
THE FULL HES/HER-PRONEURAL MODEL 

We represent Hes/Her and proneural proteins by a sin- 
gle variable in each cell. To investigate the behaviour of a 
network quantitatively, it is convenient to scale each vari- 
able such that it lies in the range [0, 1] . The full Hes/Hcr- 
Proncural model (Fig. [H^a)) can then be represented by 
the following differential equations with discrete delays: 



TABLE I: Parameter values. 
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choices of production functions that satisfy conditions 
(P and (HI). 

Numerical simulations of this model reveal a range of 
qualitatively different types of behaviour, in which oscil- 
lations can be absent, transient or sustained, and their 
phases can be locked or not (for examples, see Fig. [3]). 
To investigate the origin of these behaviours, we reduce 
the full Hes/Her-Proneural model to a variety of sim- 
pler ones (see Fig. The examination of these simpler 
networks (motifs) helps to elucidate the origins of the 
dynamics of the full Hes/Her-Proneural model. 



Th Hi = -Hi+PH{A2{t-n),Hi{t^T3)), (1) 
TaAi = -Ai+PA{Hi{t-T2),Aiit-Ti)), (2) 
Th H2 = -H2 + PH{Ai{t-Ti),H2{t~T3)), (3) 

TaA2 = -A2+PA{H2{t~T2),A2{t~n)), (4) 

where Th and Ta denote the degradation constants (the 
inverses of the linear degradation rates) of H and A, re- 
spectively . Ph and Pa are functions representing the 
rate of production of H and A respectively. The activat- 
ing and repressive action of the proneural and Hes/Her 
proteins is captured by the following constraints: 
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(5) 
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For our analysis of this and the reduced models, we 
need assume no more about the production functions 
than the conditions ^ and ([6]). For numerical simu- 
lations, specific functional forms must be assumed, and 
we take these functions to be products of increasing (/) 
and decreasing {g or h) Hill functions. Specifically, we 
assume that Pi{A,H) = fi{A)gi{H) for / = A or H, 
where 



fi{x,K,u) = — — — 
A 4- a; 

gi{x,K, v) 



(7) 
(8) 



and K and v represent the scaled threshold and the 
Hill coefficient, respectively. However, the qualitative 
behaviour of the model solutions is preserved for other 
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FIG. 3: (Color online) Examples of typical dynamics found 
in numerical simulations of the full Hes/Her-Proneural model 
shown in Fig. [2ja) and defined by Eqs. [l}{8l displaying both 
difi^erentiation and oscillations. Oscillations can be transient 
or sustained, and can be in-phase between two cells, or not. 
Values of the kinetic parameters are set the same on each 
row and listed in Table [T] Delays are given in each panel as 
(-ri,T2,r3,r4) 



IV. ANALYSIS AND SIMULATION OF 
REDUCED MODELS 

In this study, we are concerned primarily with the 
routes that cells take to differentiation. For all model 
variants, a uniform steady state exists (if* = H2, A\ = 
Aj). Biologically, this corresponds to a non-differentiated 
state (i.e. the neurogenic network is in the same state 
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in both cells). Wc therefore study the stability of this 
steady state to small perturbations, and seek to deter- 
mine the resulting dynamical behaviour of the system in 
cases where it is unstable. 

For each model variant (Fig. [2lb)-(e)), we first per- 
form linear stability analysis of the uniform steady state, 
which yields an eigenvalue equation. We then determine 
parameter values that result in the existence of neutral 
(pure imaginary) eigenvalues, which are the product of 
the imaginary unit number and the neutral angular fre- 
quency, and are associated with changes in stability. To 
confirm the nature of bifurcations associated with these 
eigenvalues, the linear analysis is followed by numeri- 
cal simulations. The eigenvalue equation for each model 
variant is derived directly from its own model equations, 
rather than by reduction of the full Hcs/Hcr-Proneural 
model. The conditions under which the full Hcs/Hcr- 
Proncural model can be reduced to the auto-activation 
two-Proneural model, and to the auto-repression two- 
Hes/Her model are discussed in Sect. IVII 

To allow systematic comparison between model vari- 
ants, we use a standard set of parameter values — which 
are listed in Table |l] — in both analysis and simulations. 
These parameter values fall into a biologically reason- 
able range [ll|, and result in representative model 
dynamics. Deviations from this standard set are noted 
explicitly. Time is measured in minutes, yielding values 
for the degradation constants that are in line with mea- 
sured values for proneural and neurogenic factors 0, [2^ . 
However, the behaviour of each model variant is also ex- 
amined with other parameter values. 

In the following analyses, 7 denotes the magnitude 
of the slope of the regulatory function {fj=i,2{x) or 
5^=1.2(2;)) at the uniform steady-state solution {H*^^ 2 
or 2)- In the cases of the Hill functions: 



7i=l,2 
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(9) 
(10) 



where, in general, K and v are different in /j=i,2(a;) and 

in gj=i^2{x), and therefore 7f4lY ^ ^(^2 ■ ^^^^ 
lowing sections, it can be seen that what determines the 
stability properties of the homogeneous steady state is 
not the precise functional forms of / and p, but the value 
of 7. 



A. Non-autonomous two-Proneural model 

The non-autonomous two-Proneural model (Fig. [2{d)) 
is represented by the following differential equations: 



TAi = -Ai+g{A2{t-T)), 
TA2 = -A2+g{Ai{t-T)). 



(11) 
(12) 



To study the stability of the uniform steady state solu- 
tion [Al = A* = A*), we set Ai{t) ^ A* + ai{t) and 



A2{t) = A* + 02 (i) in Eqs. ^ and Following lin- 

earisation, the resulting coupled equations for oi and 02 
can be uncoupled by introducing variables ai -I- 02 and 
01—02 0- The eigenvalue equations derived from the 
equations in these variables are: 

-,-At 



1 + TA = ±7 



(13) 



where the plus and minus branches are associated with 
fli — 02 and ai + 02 respectively. 

Wc consider first uniform perturbations such that ai = 
a2- In this case, only the minus branch of Eq. ()13|) is 
relevant. Assuming a pure imaginary eigenvalue A = iwc, 
with LJc real, the neutral angular frequency {lOc) is derived 
to be: 



(14) 



This eigenvalue occurs for parameter values such that 
7 > 1 and when the delay is equal to the neutral delay 



1 



= — (7r + tan-i(-rcj,)) 



(15) 



For the standard set of parameter values (Table I), the 
oscillatory period associated with the neutral eigenvalue 
(the Hopf period) — defined by Tc — 2it/uJc — is found to 
be 38.6501 min, while Tc = 13.0548 min. For r 7^ Tc, the 
minus branch of Eq. (|13p has a complex eigenvalue with 
a real part that has the same sign as t — Tc. This can be 
seen in the numerical solution of Eq. (fT3|) in Fig. ID 
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FIG. 4: (Color online) The eigenvalues (A) of the non- 
autonomous two-Proneural model (Fig. HJd)) as a function 
of delay (r). The eigenvalue equation (Eq. (|13|l ') has two 
branches, (a) On the minus branch Re{X) changes its sign 
from negative to positive as r surpasses its critical value 
(tc — 13.0548), leading to a Hopf bifurcation, (b) while on 
the plus branch it is always positive. The minus and plus 
branches therefore define respectively the oscillatory and dif- 
ferentiating properties of the system. The dotted line in (b) 
shows the approximate analytical value for the real part of 
the eigenvalue derived in Eq. (|17p . 

The linear stability analysis suggests that the uniform 
steady state becomes unstable to small uniform pertur- 
bations via a Hopf bifurcation if 7 > 1 and the delay 
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increases above the critical value. Numerical simulations 
of Eqs. pT|) and for t > are shown in Fig. [5l For 
uniform initial conditions {Ai{Q) = ^2(0)) the system 
exhibits sustained oscillations (Fig. [5ja)). This confirms 
that the neutral solution on the minus branch of Eq. (|13p 
is a Hopf bifurcation point. 
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FIG. 5: (Color online) Numerical simulations of the non- 
autonomous two-Proneural model (Fig. [IJd)). (a) With iden- 
tical initial values, Ai and A2 show sustained in-phase oscil- 
lations, (b) With slightly different initial values, Ai and A2 
show transient in-phase oscillations followed by differentia- 
tion (cf. [2^). Standard parameter values are used, as listed 
in Table 1, with r = 20 min. 

The plus branch of Eq. — which is associated with 
ai — a2 — has the same neutral angular frequency (wc) as 
the minus branch, and the neutral eigenvalue occurs for 
a neutral delay Tc- 
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was thought to be possible only in the limit of a large 
HiU coefficient {v) while in Fig. ^ v = 2. 

The positive real eigenvalue on the plus branch cor- 
responds to differentiation of the two cells (exponential 
growth of oi — 02), while the complex eigenvalue asso- 
ciated with the minus branch corresponds to oscillations 
(in fli -I- 02). Thus, if ^i(O) and ^2(0) are set slightly 
different, the behaviour of the system comprises a com- 
bination of these two fundamental modes. This can be 
seen in the numerical simulation of Eqs. (fTT|) and (fT2|) in 
Fig.lSjb), which shows a transient uniform oscillation fol- 
lowed by differentiation. Since the oscillations occur on 
the minus branch corresponding to ai -|- 02, the oscilla- 
tory profiles of Ai and A2 are in-phase. Similar transient 
oscillations leading to differentiation have been observed 
previously in a delay model of the Delta-Notch signalling 
system (Fig.IJa)) 

In general, the outcome of linear stability analysis is 
not applicable to transient system behaviour. Therefore, 
the separation of the two branches makes linear analysis 
highly valuable in this study, making possible the predic- 
tion of whether or not a transient oscillation leading to 
differentiation occurs. More striking is the fact that not 
only the initial rate of differentiation, but also the time 
taken to complete differentiation can be estimated by lin- 
ear analysis. Fig. [6] shows the time-courses of the (loga- 
rithmic) difference A2 — Ai obtained from numerical sim- 
ulations for a range of values of the delay (r) and initial 
mean values {Aq = 0.5 [^2(0) + Ai(0)]). These are com- 
pared to the growth rate predicted by the real part of the 
eigenvalue for 02—01 (i.e. the plus branch), obtained from 
linear analysis for the corresponding r (Fig. 3]), which 
provides a good estimate of the time taken to reach the 
fully differentiated state. 



Tc = — tan ^ [TlOc) , 



(16) 



which takes value 6.2702 min for the standard parameter 
set. More generally, the real part of a complex eigenvalue 
A+ satisfies 

1 + T Re{\+) = 7 e-^^(^+)^ cos(/m(A+)T). 

Since i?e(A_|_) is given by the intersection oi y ~ 
1 + Ti?e(A+) and y = 7exp(— i?e(A+)r) cos(/m(A+)r) 
on the y-i?e(A+) plane, a purely real A4- (for which 
cos(/m(A+)T) = 1) has the largest real part, and is pos- 
itive for all values of r if 7 > 1. This positive real eigen- 
value is shown in Fig. HJ^b) . By assuming a purely real 
eigenvalue (A/?+), and by applying the first-order Tay- 
lor expansion to Eq. (|13p . a simple approximate (lower 
bound) expression for is derived to be 



AH^ 



7-1 
T + -fr' 



(17) 



which is confirmed to provide a good approximation to 
the growth rate of patterns obtained in numerical simu- 
lations in Fig. IDJb) . Previously, such an approximation 



B. Auto-activation single Proneural model 

Because two sequential repressions result in a net 
activation, the non-autonomous two-Proneural model 
(Fig. [DJd)) may seem to be surrogated by the auto- 
activation single Proneural model (Fig.[2)Ie)), represented 
by the following differential equation: 



TA = -A + f{Ait~Tf)). 



(18) 



The eigenvalue equation obtained by linearisation about 
a steady state of Eq. (flSl) , is 



1 + TA = 7 e- 



-At 



which is identical to the plus branch of Eq. (|13p that rep- 
resents the differentiating nature of the non-autonomous 
two-Proneural model. It is therefore clear that the oscil- 
lation in the two-Proneural model is due to the network 
structure that two proneural proteins are mutually re- 
pressing. 
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FIG. 6: (Color online) Time taken for differentiation in 
the non- autonomous two-Proneural model (Fig. [2jd)). Time 
courses of the difference A2 — Ai — obtained from numerical 
simulations — plotted on a logarithmic scale, (a) for various 
delays (r) and initial mean values (Aq — 0.5 [^2(0) + Ai(0)]) 
with the standard parameter set, and (b) for various parame- 
ter sets with r = 1.5 Tc{H, K, v), where Tc is the critical delay. 
The time courses are compared to the growth rate predicted 
from linear analysis (dashed line), based on the real part of 
the eigenvalue (-Re(A)). 



C. Auto-activation two-Proneural model 

The auto-activation two-Proncural model (Fig. [2l^b)) 
is represented by the following dilferential equations: 

Til = -Ai+g{A2{t-Tg)) f {Alitor f)), (19) 
TA2 = -A2 + g{Ai{t-Tg)) f{A2{t^Tf)). (20) 

The eigenvalue equation for the uniform steady-state so- 
lution [Al = A*2 = A*) is derived to be: 

l+T\ = g[A*)"ffe-^^f ±J{A*)-iSe-^^^. (21) 

The purely imaginary solution X ^ iuj satisfies: 

l + B'^-C'^ + T'^uj'^ + 2B[TujsinujTf~cosL0Tf) = 0, (22) 



where B = g{A*)jf and C = f{A*)j3, 

Fig. [7] shows the neutral intercellular signalling delay 
and corresponding oscillatory period, associated with the 
pure imaginary eigenvalues, as a function of local- loop de- 
lay (ry). The eigenvalue equation (Eq. ([21]) ) is solved for 
its minus and plus branches. As tj is varied, the values of 
the neutral intercellular signalling delay and the neutral 
oscillatory period fluctuate, with the oscillatory period 
fluctuating around its value in the non-autonomous two- 
Proneural model (Fig. El^d)). A similar modulation in 
period was recently observed in a delayed coupling model 
of vertebrate segmentation (20| . Numerical simulations 
performed for the parameter values represented by the 
cross signs in Fig. [7{a) confirm that the Hopf bifurcation 
occurs on the minus branch of Eq. (PT|) . and that the 
critical intercellular signalling delay is modulated by the 
local- loop delay (data not shown). 



-minus brancli 
-plus branch 



20 40 60 80 









\ / f - 
\ / No local loop 



40 60 
local-loop delay, [min] 



FIG. 7: (Color online) Neutral intercellular signalling delay 
(a) and oscillatory period (b) associated with the pure imag- 
inary eigenvalues of the auto-activation two-Proneural model 
(Fig.HJb)), shown as a function of local-loop delay (r/). The 
eigenvalue equation (Eq. (|2ip ) is solved numerically for its 
minus and plus branches. The Hopf bifurcation occurs on 
the minus branch. Also shown in (b) is the corresponding 
Hopf period for the non-autonomous two-Proneural model, 
which lacks the local (cell-autonomous) proneural feedback 
loops (Fig. [IJd)). Numerical simulations performed for pa- 
rameter values corresponding to the cross signs in (a) con- 
firm that the Hopf bifurcation occurs on the minus branch 
of Eq. (|2H) . and that the critical intercellular signalling delay 
(Tg) is modulated by r/. 

Fig.[8]shows the results of numerical simulations of the 
model equations and (PO)) for a range of delays (r^ 
and Tf) and initial values {Ai{Q) and j42(0)). The be- 
haviour of this model variant is found to be qualitatively 
the same as that of the non-autonomous two-Proneural 
model (Fig. [21Id)): the Hopf bifurcation point exists on 
the minus branch of the eigenvalue equation (Eq. (|2ip ) 
and, since the plus branch again has a positive real 
eigenvalue, differentiation occurs when ^i(O) 7^ ^2(0). 
However, the modulatory effect of the cell-autonomous 
proneural auto-activation loop, shown in Fig. [71 is seen 
in the comparison between top and the middle right pan- 
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els. As the local-loop delay (r/) increases from to 
15, the amplitude of oscillations decreases, influenced by 
the increase of the critical signalling delay (Fig. [Tlja), 
solid curve) . It can be seen from Fig. [7| that the cell- 
autonomous proneural auto-activation loops give rise to 
'tunability' of the oscillations. 
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FIG. 8: (Color online) Numerical simulations of the auto- 
activation two-Proneural model (Fig. [IJb)). In comparison 
to Fig. [5] the behaviour is found to be qualitatively the 
same as that of the non- autonomous two-Proneural model 
(Fig. Od)). However, the modulatory effect of the cell- 
autonomous proneural loop, shown in Fig. [T] is seen in the 
comparison between top and the middle right panels. As the 
local-loop delay (r/) increases from to 15, the amplitude of 
oscillations decreases, influenced by the increase of the critical 
signalling delay (Fig. [TJa), solid curve). Values of the delays, 
and whether or not the initial values of Ai and A2 are equal, 
are specified in each panel as (rg, t/; = or 7^). 



D. Auto- repression two-Hes/Her model 

The auto-repression two-Hes/Her model (Fig. [21[c)) is 
represented by the following differential equations: 

THi = -Hi+giH2{t-Tg)) h{Hi{t-Th)), (23) 
TH2 = -H2+9{Hi{t-Tg))h{H2{t-TH)). (24) 

The eigenvalue equation for the uniform steady-state 
solution {HI = H2 = H*) has the same form as that for 
the auto- activation two-Proncural model (Eq. (PT|) ): 



(25) 



where, however, the definition of B is modified to be 
B = —g{H*)j^' because h represents a decreasing Hill 
function. 

Fig. [5] shows the neutral intercellular signalling delay 
and oscillatory period associated with the pure imagi- 
nary eigenvalues, as a function of the delay in the local 
loop. The eigenvalue equation (Eq. (|^5|) ) is solved for 
its minus and plus branches, for the first and the sec- 
ond largest periods. Unlike any eigenvalue equation of 



the three simpler models analysed so far (Fig. [S^b), (d) 
and (e)), Eq. (|25p is found not to have a neutral solution 
when T/j < 5.0483 min. 
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FIG. 9: (Color online) Neutral intercellular signalling delay 
(a) and oscillatory period (b) associated with the pure imag- 
inary eigenvalues of the auto-repression two-Hes/Her model 
(Fig. EIc)), shown as a function of the local-loop delay. The 
eigenvalue equation (Eq. (|25|) ) is solved for its minus and 
plus branches, for the first and the second largest periods 
(solid and dashed lines, respectively). This model shows the 
transition between out-of-phase and in-phase oscillations. At 
each local-loop delay, a square represents the in-phase oscil- 
lation observed in simulation with the lowest signalling delay, 
while a circle represents the out-of-phase oscillation observed 
in simulation with the highest signalling delay. The transition 
happens between them (upper panel), which are overlapping 
with a square being always above a circle. 

Numerical simulations of the model equations ([23]) and 
reveal two prominent features of the dynamics of 
this system that arc qualitatively different to those of 
the two-Proneural models (Figs. [2Kb), (d)). Fig. [TOl shows 
simulation results for a range of values of the delays (r^ 
and Th) and initial values (ffi(O) and i?2(0)). Two new 
features are apparent: first, oscillations can be sustained 
(rather than transient) even when i?i(0) 7^ -^2(0); sec- 
ond, the oscillations of Hi and H2 can be out-of-phase 
for certain values of the delays, as shown in the left bot- 
tom panel. Fig. [Tl] details the transition from out-of- 
phase oscillations to in-phase oscillations as the value of 
the critical intercellular signalling delay (r^) increases. 
Strikingly, the transition is found to be associated with 
a point-wise amplitude death [2^ . The transition found 
in numerical simulations with various is compared to 
the neutral properties estimated by the linear analysis in 
Fig. m In a cell-autonomous Hcs/Hcr oscillator, the os- 
cillatory period increases monotonically with delay p^ . 
while for the coupled cells to cycle out-of-phase, the sig- 
nalling delay must be about half a period. Therefore the 
monotonic increase in the critical Tg and in the critical 
period (Tc) suggests that the overall network behaviour 
of the auto-repression two-Hes/Her model is controlled 
by the two local autonomous Hcs/Hcr oscillators. 

In contrast, when the intercellular signalling delay (rg) 
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FIG. 10: (Color online) Numerical simulations of the auto- 
repression two-Hes/Her model (Fig. [2|c)). In comparison to 
Fig. O the behaviour is found to be qualitatively different to 
that of the non-autonomous two-Proneural model (Fig.[2jd)). 
Specifically, oscillations are sustained even with non-identical 
initial values, and oscillations can be out-of-phase as well as 
in-phase. Values of the delays, and whether or not the initial 
values of Hi and H2 are equal, are specified in each panel as 
{Tg,Th\^ or 
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FIG. 11: (Color online) Detail of the transition from out-of- 
phase to in-phase oscillations as the value of Tg increases in the 
auto-repression two-Hes/Her model (Fig. [SJc); Fig. llOf) . The 
transition is associated with amplitude death around Tg = 6. 



is varied while the local delay (r/i) is kept constant, tran- 
sitions between in-phase and out-of-phase modes occur 
sequentially, and are associated with modulations in am- 
plitude and period (Fig. I12p . The mode transitions re- 
peat as Tg is increased, but not in an identical manner. 
Unlike the first transition shown in Fig. 1111 a clear am- 
plitude death is not observed at higher-order transitions. 
Further examination reveals that the sequential transi- 
tions are not induced by the network structure of two 
mutually repressive autonomous oscillators, but by the 
structure of two mutually influencing autonomous oscil- 
lators, forming an overall positive feedback loop. In- 



deed, when intercellular interactions are modelled by an 
increasing Hill function that represents activation, the 
transition still occurs, but in a reverse order, starting 
from the in-phase mode at Tg = (data not shown). 
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FIG. 12: (Color online) Dependence of the amplitude (top 
panel) and period (middle panel) of oscillatory solutions on 
the intercellular signalling delay (rg) in the auto-repression 
two-Hes/Her model. The amplitude of an oscillatory solu- 
tion is defined as the half difference between peak and trough 
values. The observed modulations of amplitude and period 
reflect the transition between out-of-phase (OP) and in-phase 
(IP) oscillation modes; the comparison between lag (time dif- 
ference between oscillation peaks in neighbouring cells) and 
the half-period of oscillations shown in the middle panel illus- 
trates the sequential switching between out-of-phase and in- 
phase oscillation modes. The bottom four panels show time 
courses obtained by numerical simulation of Eqs. (|23|) and 
(|24p at the four values of Tg marked in the top panel. The 
first transition corresponds to Fig. 1111 The dashed lines in 
the top and the middle panels represent respectively the am- 
plitude and the period in the single-cell Hes/Her model with 
local delay — 15. 



The intercellular coupling is found to have additional 
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dynamical consequences. While an isolated cell contain- 
ing the Hes/Her feedback loop cannot exhibit sustained 
oscillations when the delay (r/j) is below its critical value, 
such a local- loop critical delay (t/j cr) can be lowered by 
the intercellular coupling, as shown in Fig. I13i where 
T/i = 10 < Th cr = 13.0548. Furthermore, oscillation 
death is found to happen in a definite range of the inter- 
cellular signalling delay (vg). 
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FIG. 13: (Color online) Coupled sub-critical oscillators illus- 
trating a finite- Tg-range oscillation death. The inset in the 
top panel shows the oscillation death in detail. 



relations between the neutral delay in the local loop 
(Tfi cr) and the intercellular signalling delay (xg) shown 
in Fig. [T31 for both the out-of-phase (green) and the in- 
phase (blue) modes. For example, for t/j = 10 (Fig. [T5)l . 
when < Tg < 2 only the out-of-phase mode is al- 
lowed, whereas when Tg > 5.1 only the in-phasc mode 
is allowed. In between {2 < Tg < 5.1), no oscillation 
is allowed, explaining the origin of the oscillation death 
shown in Fig. [131 However, the comparison of periods 
to simulation results in the lower panel shows that as 
Tg increases, more and more branches appear, and the 
model behaviour becomes more non-linear. The actual 
transitions observed in simulations do not coincide with 
the boundaries between the regions that linear stability 
analysis predicts to be competent (solid) and incompe- 
tent (dashed) for sustained oscillations. 
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The three features observed in numerical simulations — 
the oscillation in an originally subcritical delay range 
(Th < Th cr), the transition between out-of-phase and 
in-phase oscillatory modes, and the finite- r^-range oscil- 
lation death — can be understood analytically as follows. 
Linearisation of Eqs. ([23]) and ([24]) around the uniform 
steady state yields: 



20 




(26) 

where a = g{H*)j''e-^'''\ (3 = h{H*)-i^e~^^^ . The asso- 
ciated eigenvalue equation, given by dctM = 0, is: 

l + T\ = -a±l3, (27) 



which is the same as Eq. ([^5]) . By substituting Eq. ([?7|) 
into Eq. ([^5]) , it can be seen that the plus branch is associ- 
ated with the out-of-phase oscillation mode: {Hi,H2) = 
(1,-1), whereas the minus branch is associated with 
the in-phase oscillation mode: {Hi,H2) = (lil)- Nu- 
merical evaluation of the eigenvalue equation yields the 
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FIG. 14: (Color online) The origin of oscillation death, and 
the transition between the out-of-phase and the in-phase os- 
cillatory modes in the auto- repression two-Hes/Her model 
(Fig. [3c)). Neutral local- loop delay (upper panel) and 
oscillatory period (lower panel) are shown as functions of 
the intercellular signalling delay Tg. The eigenvalue equa- 
tion (Eq. ([25}) is solved numerically for its minus and plus 
branches. Neither the out-of-phase branch (green) nor the 
in-phase branch (blue) have a solution for tu = 10 for 
2 < Tg < 5.1, corresponding to the oscillation death observed 
in Fig. 1131 Simulated periods for Th = 15 (shown also in 
Fig. 1121 middle panel) are shown by blue and green dots in 
the lower panel, where the neutral period obtained by lin- 
ear stability analysis is shown by solid and dashed curves 
(correspond to the neutral delay being smaller or larger than 
Th = 15, respectively). Note the different horizontal scales in 
the two panels. 



It is striking that only by changing the local loop at- 
tached to the non-autonomous two-X model from auto- 
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activation (Fig. [2jb) where X is Proneural) to auto- 
repression (Fig. [2Kc) where X is Hes/Her), the network 
behaviour changes completely. This simple example, mo- 
tivated by the structure of the neural differentiation net- 
work, highlights the importance of examining in detail 
the structure of any gene regulatory network, even if it 
is centred on a small and seemingly simple motif. 

Richness in the behaviour of this auto-rcprcssion two- 
Hes/Her model can be seen in Fig. [151 where the largest 
neutral period (Tc) associated with the eigenvalues is 
plotted for different values of thee threshold {Kf^) and 
Hill coefficient {vh) in the local feedback loop. As Kh 
increases, the behaviour of this model becomes simi- 
lar to that of the auto-activation two-Proneural model 
(Fig. E^b)), while as i^h decreases, the behavour becomes 
similar to that of the non- autonomous two-Proneural 
model (Fig.IKd)). 
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EXTENSION TO iV-CELL RING MODEL 



The two-cell models discussed so far can be naturally 
extended to one-dimensional arrays of cells. To avoid po- 
tential boundary effects, we consider the specific case of a 
ring of N cells labelled with a single index i = 1, 2, . . . , 
(i.e. a line of cells with periodic boundary conditions im- 
posed). As an example, we study the auto-repression 
Hes/Her model. We assume that each cell signals equally 
to both its nearest neighbours, yielding the following 
model equations: 



(28) 

where i denotes cell number and the imposition of peri- 
odic boundary conditions implies that iJ-i = Hn and 
i^AT+i = Hi. Linearising Eq. (j28p around the uniform 
steady-state solution {Hi ~ H*) by expanding Hi as 
Hi = H* + HiC^* yields the following eigenvalue equa- 
tions: 



= AH.^i + {B- XT)H, + AH,+i, i = l,2,...,N, 

(29) 

where 



FIG. 15: (Color online) The largest neutral period (Tc) associ- 
ated with the eigenvalues of the auto-repression two-Hes/Her 
model (Fig. [2jc)) for different values of the threshold {Kh) 
and the Hill coefficient {uh) in the local feedback loop. As 
Kh increases, the behaviour of this model becomes similar to 
that of the auto- activation two-Proneural model (Fig. [2jb)), 
while as Uh decreases, the behavour becomes similar to that of 
the non- autonomous two-Proneural model (Fig. [2jd)). Note 
the different vertical scales in each plot. The dashed-dot 
line in each panel represents the critical period in the non- 
autonomous two-cell model. Varied parameters are given in 
each panel as {Kh, Vh)- 



A{\) = --h{H*)^Se-^^^, 
B{X} = -I- g{H*)-f^e~^''K 



Eq. can be represented in matrix form as: 



/0\ I B -XT A 00 

A B-\T A 00 

A B-\T A Q 



\q J \ A 00 









A 





H2 

H3 



A B-XTj \h^J 



( Hi \ 

H2 
H3 



\Hn J 
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where it is important to note that A and B arc non- 
polynomial functions of A. The eigenvalues A are deter- 
mined by det M = 0. 

We first note that the phase relations between adjacent 
cells in oscillatory solutions can be determined from the 
form of the eigenvector {Hi, H2, ■ ■ ■ Hn) corresponding 
to each eigenvalue. For any value of N, Eq. ([50]) has an 
eigenvector (1, 1, . . . , 1) with an eigenvector determined 
by the solutions of TA = B + 2A. This corresponds to 
an oscillatory mode where all cells are in-phase. Further- 
more, for any even value of iV, Eq. pO|) has an eigenvec- 
tor (1,-1,1,. ..,1,-1) with an eigenvector determined 
by the solutions of TA = B — 2A. This corresponds to an 
oscillatory mode where all cells are out-of-phase. These 
two cases are simple extensions of the dynamics observed 
in the two cell model. 

To illustrate potential extensions to the dynamics ob- 
served in the two cell model, we consider the specific case 
iV = 4: 



f B- XT A A \ 

A B - XT A 
A B- XT A 

\ A A B -XT J 



(31) 



and 



detl 



{-2A + TX- B)i2A + TX- B){-B + TXf 



Therefore, the eigenvalues are determined by the solu- 
tions of: 



peaks per oscillatory period and "intermittent" oscilla- 
tions being common. An example for A'^ = 7 is shown 
in Fig. [iW b). These features are also observed in simu- 
lations on regular square and hexagonal two-dimensional 
arrays of cells (results not shown). 
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The explicit expression of Eqs. ([5^ and are found 
to be the same as the eigenvalue equations derived for 
the auto-repression two-Hes/Her model (Eq. (|25|) ). and 
for the auto-repression single Hes/Her model (the minus 
branch of Eq. p^ ). respectively. 

As above, the phase relations between adjacent cells in 
oscillatory solutions can be determined from the form of 
the eigenvector (Hi, H2, H-s, H4) corresponding to each 
eigenvalue. From Eq. (|3ip . the eigenvector on the minus 
branch of Eq. ((32l) is found to be (1, —1, 1, —1), represent- 
ing the out-of-phase oscillatory mode, while the eigen- 
vector on the plus branch is (1,1,1,1), representing the 
in-phase mode. In contrast, the eigenvector associated 
with Eq. ((33|) . which is present when A^ = 4, 8, 16, . . ., 
is found to be (1, 0, —1, 0), representing amplitude death 
for every other cell in the array, with the remaining cells 
being out-of-phase. This is a new dynamical feature that 
the two-cell model fails to capture. An example of this 
mode for A^ = 4 is shown in Fig. [TOT a). 

For larger values of A^, interaction between multiple 
modes can result in more complex oscillatory behaviour, 
in which waves of phase and amplitude differences can 
propagate around the ring. Furthermore, the oscilla- 
tory profile for each cell is often complex, with multiple 



FIG. 16: Complex oscillatory dynamics in one-dimensional 
arrays of coupled cells. Simulations of Eqs. (|28|l , representing 
the the auto-repression Hes/Her model in a ring of A'' cells. 
Simulations were carried out using the standard parameter 
values (Table 1) with — 20, Tg = 8 and with initial his- 
tory Hi{t) — 0.1 + O.l^i for —msLx{Tg,Th) < i < 0, where 
the are drawn from independent uniform distributions on 
[0, 1] (a) Amplitude death combined with out-of-phase oscil- 
lations: N = 4. (b) Complex dynamics resulting from mode 
interactions: N = 7. 



VI. ANALYSIS AND REDUCTIONS OF THE 
FULL HES/HER-PRONEURAL MODEL 

In this section, from the viewpoint of eigenvalue equa- 
tions, it is discussed how the full Hes/Her-Proneural 
model (Fig. HJa)) is related to the auto-activation 
Proneural model (Fig. [2Ub)), and to the auto-repression 
Hes/Her model (Fig. [Sic)). 

For the uniform steady-state solution {H'^ = = H* , 
Al = A2 — A*) of the full Hes/Her-Proneural model, the 
eigenvalue equation is derived from its model equations 
©-dll) to be: 



{l+THX+ae-^^'){l+TAX 



- At4 ^ 



±Me"^^*°S (34) 
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where a = /l(A*)7^^ b = /2(A*)7^'S c = gi{H*)^f\ 
d = 92{H*)-ff^ , and not = n + ra. 

First, explicit forms of the neutral angular frequency 
(lUc) and delay {rtot) are derived for a simple example, 
where Ta = Th = T, T3 = T4 = Ta, and = /a = 
gi = §2- The last assumption means that K and v are 
respectively the same all in /i, /2, gi, and (72- Then, 
from Eqs (9) and (10), A* = H* , and consequently a = 
c = b = d -.^ (f). Therefore, Eq. (134]) becomes: 

For a neutral solution, Eq. ([55]) becomes: 

(1 - T^cj^ - c/)^ cos 2ujTaf + {2Tuj + 0^ sin2wTa)^ = 0'', 

(36) 

tanwTtot = TpT^ -73. 7, TA. (37) 

1 - T-^UJ'^ - 0^ cos 2ujTa 

For the plus branch in Eq. ([M|) : 

Tfot = - max[ tan"i(-A), tt + tan"i(-A) ], (38) 
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FIG. 17: (Color online) Neutral (a) intercellular signalling 
delay and (b) period, associated with the pure imaginary 
eigenvalues of the full Hes/Her-Proneural model, shown as 
a function of the local-loop delay. The eigenvalue equation 
(Eq. psp ) is solved for its minus and plus branches, for the 
first and the second largest periods (solid and dashed lines, 
respectively). In comparison to Figs. [7] and [O] the auto- 
repression two-Hes/Her circuit component is found to be dom- 
inant in this full Hes/Her-Proneural model. 



whereas for the minus branch in Eq. (|34p : 

Ttot = - max[ tan"^(+A), tt + tan"^(+A) ]. (39) 

LO 

Fig. [17] shows the neutral (a) intercellular signalling 
delay and (b) period, as a function of the local-loop delay 
{to)- The eigenvalue equation (Eq. ([35])) is solved for 
its minus and plus branches, for the first and the second 
largest periods. In comparison to Figs. [7] and [51 the auto- 
repression two-Hes/Her circuit component is found to be 
dominant in this full Hes/Her-Proneural model, although 
the neutral values are found to exist for any Ta , unlike in 
the auto-repression two-Hes/Her model, where no pure 
imaginary solution can exist for Tf less than 5.0483 [tj 
corresponds to in this section). 

In the following, based on the eigenvalue equation 
of the full Hes/Her-Proneural model (Eq. ([M])), we 
clarify the conditions under which this full Hes/Her- 
Proneural model can be reduced to the auto-activation 
two-Proneural model (Fig. [Hb)), and to the auto- 
repression two-Hes/Her model (Fig. [2l^c)). 



A. Reduction to the auto-activation two-Proneural 
model 

The eigenvalue equation of the full Hcs/Hcr-Proneural 
model (Eq. ([M)) ) is reduced to the eigenvalue equation 
of the auto-activation two-Proneural model (Eq. ([H])) 
when: [a] Th = 0; [b] = 0; [c] H* = A*- [d] 
g2{H*)jf'^ ~ 1; and [c] ri = 0. These conditions must 
hold irrespective of whether fi and gi are Hill functions 
or not. If Hill functions f{x, K, v) and g{x, K, v), where 



K and v are the scaled threshold and the Hill coeffi- 
cient, are employed, condition [b] can be met by set- 
ting K in 52 to be +00. Because this operation makes 
g2{H*) = 1, from the steady-state solution of the full 
model {H* = ./i(A*) g2[H*)), condition [c] is found to 
require that /i is not a Hill function but a linear func- 
tion: fi{x) = X, whereby 7-^^ = 1, making all five con- 
ditions satisfied. The last operation {fi{x) = x) means 
that for a sequence of two functions / and g to be rep- 
resented only hy g, f needs to be f{x) = x. In short, 
the auto-activation two-Proneural model (Eqs. ([TO]) and 
([20)) ) is obtained by assuming the following conditions 
on the full Hes/Her-Proneural model (Eqs. ([Il)-([1|)): the 
activation from a proneural protein (^^=1^2) to the adja- 
cent Hes/Her (i?i=2,i) is linear and instantaneous; and 
Hes/Her is not associated with an auto-repression loop. 



B. Reduction to the auto- repression two-Hes/Her 
model 

The eigenvalue equation of the full Hes/Her-Proneural 
model (Eq. ([M)) ) is reduced to the eigenvalue equation 
of the auto- repression two-Hes/Her model (Eq. ([^S]) ) 
when: [a] Ta = 0; [b] 7/= = 0; [c] h{A*) = gi{H*)- 
[A] f2{A*)"ff^\ [e] Ti — 0. These conditions must hold 
irrespective of whether fi and gi are Hill functions or 
not. If Hill functions are employed, condition [b] can be 
met by setting K in /2 to be 0. Because this operation 
makes f2{A*) ~ 1, from the steady-state solution of the 
full model [A* ~ gi{H*) /2(A*)), condition [c] is found to 
require that /i is not a Hill function but a linear function: 
fi{x) = x, whereby 7-^^ = 1, making all five conditions. 
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VII. GENERAL DISCUSSION 

The development of multicellular organisms is a pro- 
cess of sequential and concurrent cell differentiations, the 
timings of which must be tightly controlled. Recent ex- 
perimental studies have revealed that in the case of verte- 
brate neural differentiation, cell differentiation can occur 
after a transient oscillation in cell state [23|- Since neu- 
ral differentiation in vertebrates occurs over a consider- 
able time interval, these transient oscillations may play a 
role in the scheduling of neural differentiation in relation 
to other developmental events. The occurrence of such 
transient oscillations on the route to neural differentia- 
tion was previously predicted in a simple delay model of 
Delta- Notch- medated lateral inhibition p6j . 

In the present study, we have analysed a more detailed 
model of the neural differentiation network. The model 
has a nested-loop structure, with intercellular signalling 
(as mediated by interactions between DSL ligands and 
Notch family receptors) coupled to local cell-autonomous 
feedback loops. Using a combination of stability analy- 
sis and numerical simulation, we have shown that the 
incorporation of local feedback loops has potentially sig- 
nificant impacts on the dynamics of the signalling net- 



work. In particular, the time delays in the local feedback 
loops were found to play a central role in controlling the 
behaviour of the whole network: whether it heads to- 
wards differentiation; whether it shows an oscillation; or 
whether such an oscillation is sustained or transient, as 
well as providing tunability in the amplitude and period 
of oscillations. 

Many classes of auto-regulatory genes have been identi- 
fied. Because these genes are also nodes in larger genetic 
regulatory networks, nested feedback loop structures are 
rather common. The results of the present study high- 
light the need for careful examination of the predictions 
made by non-delay network models, which include the 
Drosophila neurogenic model studied by Meir et al. that 
yielded a conclusion that the total system was robust to 
local changes to the network circuitry [l3| . 

In a more general perspective, a biological system is 
known to involve a relatively small number of genes, hav- 
ing numerous features and functions. Some of the new 
functions may have been acquired by the addition of a 
few local loops to the old conserved ones. Moreover, be- 
cause the loops considered in the present study are all 
very simple, the outcomes of the present study may have 
relevance to non-biological systems as well. 
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